
#include<stdiostream.h>
#include<string.h>
#include<stdlib.h>
#include<ctype.h>



#include <vector>
#include<complex>
using namespace std;



//begin Power spectrum code

class PowerSpectrum  
{
public:
	PowerSpectrum();
	virtual ~PowerSpectrum();
	char * strCRegion;
	double *s;
	double *x;
	int  doPowerSpectrumAnalysis(int); 
	void	setSequence(const char * seq);
	
};


PowerSpectrum::PowerSpectrum()
{
	s=NULL;
	x=NULL;
	strCRegion=NULL;

}

PowerSpectrum::~PowerSpectrum()
{
	
	delete []s;
	delete []x;
    delete []strCRegion;

}

void PowerSpectrum::setSequence(const char * seq)
{

delete []strCRegion;

strCRegion=new char[ strlen(seq)+1 ];

memset(strCRegion, '\0' , strlen(seq)+1 ) ;

strcpy(strCRegion,seq);

}


int PowerSpectrum::doPowerSpectrumAnalysis(int pc) 
{


const double PI = 3.14159274;

complex<double>  i1(0.0,1.0);

long lenStr = strlen(strCRegion);


delete []s;
delete []x;


//cout<<"Coding Region Length:"<<lenStr<<endl;


s	= new double[lenStr/2 + 1];
x	= new double[lenStr/2 + 1];

char nuclt[5]= "ATGC";
int const lenWindow = (lenStr/pc)*pc;
 
double li=(1.0/lenWindow);

double winsq		= li+li*li;
//long  wincube	    = lenWindow*lenWindow*lenWindow;


int countA =0;
int countC =0;
int countG =0;
int countT= 0;

complex<double> s1(0.0,0.0);

double sq=0.0;
double sq1=0.0;
double s2=0.0;
double sbar =0.0;
double p=0.0;

	
//ofstream out("power");

	for(long i =0; i < lenWindow; i++ )
	{

	if( strCRegion[i] == 'A') 
		{ countA++; continue;}
  	if( strCRegion[i] == 'G') 
		{ countG++; continue;}
  	if( strCRegion[i] == 'C') 
		{ countC++; continue;}
  	if( strCRegion[i] == 'T') 
		{ countT++; continue;}

	} 

//out<<"Count A  :"<<countA <<"    countT :"<<countT <<"countG   :"<<countG <<"countC   :"<<countC <<endl;

 
sq1 = 0.0+ pow(countA,2)+pow(countG,2) + pow(countT,2) + pow(countC,2);

double f =0.0;



long k;
for(k = 1 ; k <=lenWindow/2  ; k++ )
{

	s2	=0.0;
	sbar	=0.0;
	p	=0.0;
	f= ((1.0)*k)/lenWindow;

// for each Nucleotide. A T G C


	for(int i = 0; i <4 ; i ++ )
	{
		s1= complex<double>(0.0,0.0);


		for (int bi=0; bi<lenWindow; bi++)
		{
				if(strCRegion[bi] == nuclt[i] )
					s1= s1+exp(-i1*2.0*PI*f*(bi+1.0));
		}

		s2 = s2+ (abs(s1)*abs(s1));


	}
	
//	to calculater power specturm and code region 

		sq   = sq1*li*li*li;
	     
		sbar =(winsq-sq);
		sq   = sbar * (lenWindow/2.0);

// to calcuater signal to noise ratio

		s2 = s2/(sq*lenWindow*lenWindow);
		
		p= (s2/sbar);
		
		
	s[k-1]=f;
	x[k-1]=p;
	

//	out<<"\n Frequency:"<<f <<" Power  : "<<p;

} // outer for loop



int numPeaks=0;
int genePeak=0;
//WindowAnalysis wobj; 

for( k = 1 ; k <= (lenWindow/2)  ; k++ )
{

	if( x[k-1] > 4.0)
	{
	   //if(  s[k-1] == 1.0/3.0  )] //added by Sreenu on mandars suggestion on 19/04/2002
	  //if( ( s[k-1]>= 0.32 ) && ( s[k-1] <= 0.34 ))
             
		if(s[k-1]==(1.0/pc)) 
		  genePeak=1;
	        else
	         numPeaks++;
	}

}


//delete []s	  ;
//delete []x	  ;
//s =NULL  ;
//x =NULL  ;


if(genePeak==1)
{  
         if(numPeaks==0)
	    return 1;
	 else
	   return 2;
}
else
{
        return 0;
}



}  // end of member function.

//end Power spectrum code


//begin Windowanalysis code

struct structCR
{
	long start;
	long end;
	char type;

};


class WindowAnalysis //__declspec(dllexport)   
{
public:

		
	long    *s;
	double   *x; 
    int     pcity; 
	
	vector<structCR> vecCRegion;
			 WindowAnalysis();
			 WindowAnalysis(int len ,double thold,long intron , long exon);
	virtual ~WindowAnalysis();
	void	  setWindowLen(int len);
	void	  settherHold(double thold);
	void      setIntronExon(long intron , long exon );
	int		  doFt(const char  * strSeq , char *logfile);
	void      setpcity(int pc);


private :

	int		lenWindow;
	double  thersHold;
	long	exonCutOff;
	long	intronCutOff;
	
	int		findexp(int WindowSize);
	long    framecount;
	complex<double>  * expval;
	WindowAnalysis& operator =(const WindowAnalysis &cwa);
	WindowAnalysis(const WindowAnalysis &cwa);

};

WindowAnalysis::WindowAnalysis()
{
	expval = NULL;
        s      = NULL;
        x     = NULL;
        framecount=0;


}

WindowAnalysis::~WindowAnalysis()
{
	
	delete []expval;	
	delete []s	  ;
	delete []x	  ;
	printf("Destructor of WindowAnalysis ");
	fflush(stdout);
  //delete []CRegion;
	//this->vecCRegion.clear();


}

WindowAnalysis::WindowAnalysis(int len , double thold ,long intron , long exon)
{

	lenWindow=len;
	thersHold=thold;
	exonCutOff=exon;
	intronCutOff =intron;
	
	s	   = NULL;
    x	   = NULL;

	framecount=0;
//	cout<<"ThresHold Value "<<thersHold;
//	cout<<"Intron :"<<intronCutOff<<" Exon :"<<exonCutOff<<endl;

    expval=new complex<double>[lenWindow];

}


void WindowAnalysis::settherHold(double thold)
{
	thersHold =thold;
//	 cout<< "ThresHold Value "<<thersHold;

}
void WindowAnalysis::setpcity(int pc)
{
	pcity=pc;
}

void WindowAnalysis::setWindowLen(int len)
{
    lenWindow=len;
	delete [] expval;
	expval=new complex<double>[lenWindow];

}
void  WindowAnalysis::setIntronExon(long intron , long exon )
{
	 exonCutOff=exon;
	 intronCutOff =intron;


}
int WindowAnalysis::doFt(const char  * strSeq , char * logfile)
{


const double F3 = 1.000/pcity;
const double PI = 3.14159274;


//ofstream	out("output" );
//ofstream	out(logfile ,ios::app);

long lenStr = strlen(strSeq);

if(lenStr<lenWindow)
{
   //out<<"Sequence length is less than windowlength."<<endl;
   return 1;
}

//cout<<"\n Sequence Length :"<<lenStr<<endl;

int result= findexp(lenWindow);

//long *CRegion = new long[lenStr-lenWindow+ 1];
//long *s	    = new long[lenStr-lenWindow+ 1];
//double *x	= new double[lenStr-lenWindow+1];

//delete [] s;
//delete [] x;
s	    = new long[lenStr-lenWindow+ 1];
x	= new double[lenStr-lenWindow+1];


char nuclt[5]= "ATGC";

 double li=(1.0/lenWindow);

 
 double winsq		= li+ (li*li);

double winsq1		= (1.0/lenWindow)+(1.0/(lenWindow*lenWindow));
double wincube	    = lenWindow*lenWindow*lenWindow;


int countA =0;
int countC =0;
int countG =0;
int countT= 0;

complex<double> s1(0.0,0.0);

double sq=0.0;
double sq1=0.0;
double s2=0.0;
double sbar =0.0;
double p=0.0;

framecount=0;



for(long k = 0 ; k < lenStr-lenWindow+1 ; k++ )

{
	sq		=0.0;
	s2		=0.0;
	sbar	=0.0;
	p		=0.0;

	// Frame wise  count of A C G T seperately


if(k==0)
{
	for(long i =k; i < k+lenWindow; i++ )
	{

	if( strSeq[i] == 'A') 
		{ countA++; continue;}
  	if( strSeq[i] == 'G') 
		{ countG++; continue;}
  	if( strSeq[i] == 'C') 
		{ countC++; continue;}
  	if( strSeq[i] == 'T') 
		{ countT++; continue;}

	} //end of count for look.

}else
{

switch(strSeq[k-1])
{
case 'A' :
    { countA--; break;	}
case 'T' :
    { countT--; break;	}
case 'G' :
    { countG--; break;	}
case 'C' :
    { countC--; break;	}
}

switch(strSeq[k+lenWindow-1])
{
case 'A' :
    { countA++; break;	}
case 'T' :
    { countT++; break;	}
case 'G' :
    { countG++; break;	}
case 'C' :
    { countC++; break;	}
}

 
}

//out<<"Count A  :"<<countA <<" Frame count :"<<k<<endl;
//out<<"Count T  :"<<countT <<" Frame count :"<<k<<endl;
//out<<"Count G  :"<<countG <<" Frame count :"<<k<<endl;
//out<<"Count C  :"<<countC <<" Frame count :"<<k<<endl;

sq1 = 0.0+ pow(countA,2)+pow(countG,2) + pow(countT,2) + pow(countC,2);

// for each Nucleotide. A T G C


	for(int i = 0; i <4 ; i ++ )
	{
		s1= complex<double>(0.0,0.0);


		for (int bi=0; bi<lenWindow; bi++)
		{
				if(strSeq[k+bi] == nuclt[i] )
		 			s1= s1+expval[bi];
		}
		//cout<<s1<<endl;

		s2 = s2+ (abs(s1)*abs(s1));

	}


//	to calculater power peak 

		
	    sq   = sq1*li*li*li;
		sbar =(winsq-sq);
		sq   = sbar * (lenWindow/2.0);

// to calcuater signal to noise ratio

		s2 = s2/(sq*lenWindow*lenWindow);

		p= s2/sbar;


	s[framecount]=k+(lenWindow/2);
	x[framecount]=p;

	//out<<s[framecount]<<"                "<<x[framecount]<<"\n";
    framecount++;


} // End of Outer for loop(Frame wise)

//return(0);
//} 
// end of Ft function..
//void WindowAnalysis::findCodingRegion(long seqfileId , const char *  seqId)
//{

long *CRegion = new long[framecount+1];
int    flag=0; 

long i =0 , m =0;

	for (m =0; m<framecount; m++)
	{
       
		if(x[m]>=thersHold)
		{
			if(flag == 0)
			{
				
				CRegion[i] =s[m];
				//out<< "\n START :" <<s[m]<<" : "<<x[m];
				//i++;
				flag =1;
			}

		}

        else
        {
            if(flag==1)
       	     {
               if(x[m]<thersHold)
               {
				 	
					//out<< "\n START :" <<CRegion[i];
					//out<<"\t\tEND: "<<s[m-1]<<" : "<<x[m-1];		 
					i++;
					CRegion[i]=s[m-1];
					i++;
				}
				flag=0;

			 }
			
        } // end of if.
	}// End Coding region loop. 



if (flag==1)
{
	//out<< "\n START :" <<CRegion[i];
	//out<<"\t\tEND: "<<s[m-1]<<" : "<<x[m-1];		 
	i++;
	CRegion[i]=s[m-1];
	i++;
}

//delete []CRegion;
// that data is deleted 07/02/02 To inprove performance.

//delete []s	  ;
//delete []x	  ;
//s=NULL;
//x=NULL;
 
//out<<"\n Final Coding Regions   :"<<endl;

flag =0;
int code =0;
long sta=0 , stp =0;
structCR obj;
obj.type=' ';
					
for(int j=0 ; j<i ; j +=2)
   {
	if(flag ==0)
	{
		sta = CRegion[j];
		stp = CRegion[j+1];
        if ( (stp-sta)>= exonCutOff ) code =1; else code =0; 	
		flag = 1 ;
    }
	else
	{                      
		if( (CRegion[j]-CRegion[j-1]) >=intronCutOff  )
		{
			if(code ==1)
			{	
				//out <<"\nStart :"<<sta <<"  End  :" << stp;  			
				// inset in database  start new set.
				    obj.start=sta;
					obj.end=stp;
					vecCRegion.push_back(obj);
				//insetcodes( sta , stp,seqfileId,seqId);
				sta = CRegion[j];
				stp = CRegion[j+1];
			
				if ( (stp-sta)>= exonCutOff ) code =1; else code =0; 	
				flag = 1 ;
			}
			else
			{ 
				//cout << forget the old set.
				  if ((stp-sta)>= exonCutOff ) 
				  { 
					//out <<"\nStart :"<<sta <<"  End  :" << stp;
					obj.start=sta;
					obj.end=stp;
					vecCRegion.push_back(obj);
					//insetcodes( sta , stp,seqfileId,seqId);
				  }
				// may need to be  commented later.
				sta = CRegion[j];
				stp = CRegion[j+1];
				if ( (stp-sta)>= exonCutOff ) code =1; else code =0; 	
				flag = 1 ;
			}
		}
		else
		{

        	if( ( CRegion[j+1]-CRegion[j] ) >= exonCutOff ) code =1;
			stp = CRegion[j+1];
        	flag =1;

		}
	}
 }

if ((code == 1) || ((stp-sta) >exonCutOff) )
{

	//out <<"\nStart :"<<sta <<"  End  :" << stp;
	obj.start=sta;
	obj.end=stp;
	vecCRegion.push_back(obj);
	//mapCRegion[sta] =stp;
 	//insetcodes( sta , stp,seqfileId,seqId);		
}


//delete []CRegion; 
//CRegion=NULL;


vector<structCR>::iterator it;

PowerSpectrum powerObj;

char * strtemp=NULL;



 for(it =vecCRegion.begin(); it !=vecCRegion.end(); it++)
 {
	//long start =it->start;
	//long end =it->end;

	strtemp = new char[( (it->end - it->start) +1)];

	memset(strtemp , '\0' ,(it->end - it->start) +1 );

	strncpy( strtemp , &strSeq[it->start] ,(it->end - it->start) );

	//strncpy( strtemp , &strSeq[start] , end -start);

    //printf("Length : %d \n ", strlen(strtemp));

	//printf("start %ld , End  %d  length %d\n " , it->start ,it->end , (it->end - it->start));
	//fflush(stdout);

	powerObj.setSequence(strtemp);
	int peaks =powerObj.doPowerSpectrumAnalysis(pcity);
	if(peaks==0)
	{
		it=vecCRegion.erase(it);
		if( it !=vecCRegion.end() )
		--it;
		else
			break;
	}
	else if(peaks==1)
	{
		it->type ='S';
	}
	else if(peaks==2)
	{
		it->type ='M';
	}
 delete []strtemp;
 strtemp = NULL;
 }
 


  //store the coding regions in an array


//int elsize;
 //elsize=vecCRegion.size();



delete []strtemp;
strtemp = NULL;


//delete []CRegion; 
//CRegion=NULL;

//delete []s	  ;
//delete []x	  ;

return(0); 

} // end of Coding region finding.
 

// to find exponential positional values to window size .

int WindowAnalysis::findexp(int WindowSize)
{


complex<double>  i1(0.0,1.0);
const double PI = 3.14159274;
const double F3 = 1.000/3;

for(int k=0; k<WindowSize; k++)
{ 
    	
	expval[k] = exp(-i1*2.0*PI*F3*(k+1.0));		
	
	//cout<<cos(-1*2.0*PI*F3*(k+1.0))<<endl;
}

return(0);

}



// end Window Analysis code 





class ftcds
{

	
	WindowAnalysis  *winobj, *winobj1;
	PowerSpectrum pobj;


	FILE *fp1,*fp2,*fp3,*fpl,*fpc;
	char *seq,*seqs;//store the sequence to find the PNR;
	long begpos,endpos;
	//ifstream begptr,endptr;
 
   

public:

	void doGenscan(char *,char*,int,int,int,int,int);

	ftcds()
	{
	   //seq = NULL;
	   winobj = new WindowAnalysis();
	   winobj1 = new WindowAnalysis();
	   
	
	}

	~ftcds()
	{
		delete [] seq;
		seq = NULL;
		cout<<"Destructor of ftcds"<<endl;
		delete winobj;
		delete winobj1;
		
	}


};

 

 void ftcds::doGenscan(char *inputfile,char *outfile,int winlen,int thold,int pcity,int exoncut,int intcut)
 {
	 
	 //variable where the sequence will be stored

     int seqno=1;
	 	 
	 //char fname[50];
	 char ch,ch1;
	 char *seqname;
	 int snin;
     
	 vector<structCR>::iterator it1;
     
	 
	 fp1=fopen(inputfile,"r");

	 if(fp1==NULL)
	 {
		cout<<"File doesn't exist"<<endl;
		exit(1);
	 }

	 //create a temporary file to store the sequence without the newline characters
	 	 
	 fpc=fopen(outfile,"w");
	 if(fpc==NULL)
	 {
		cout<<"Could not open output file"<<endl;
		exit(1);
	 }
	 
	 ch=fgetc(fp1);
	 
	while(ch!=EOF)
	{
		
	   	   
		fp2=fopen("tmpseq.tmp","w+");
	 
		if(fp2==NULL)
		{
			cout<<"Could not create temporary file"<<endl;
			exit(1);
		}
						 
		do 
		{
		 
		 //check for the initial '>'symbol put the characters after the next line

			if(ch=='>')
			{

				begpos=ftell(fp1);

				 while(ch!='\n')  
				 {
		 			 ch=fgetc(fp1);
				 } 

				 endpos=ftell(fp1);
				 seqname=new char[endpos -begpos + 2];
				 memset(seqname,'\0',endpos-begpos  +2);

				 fseek(fp1,begpos,SEEK_SET);
				 snin=0;
				while(true)
				{
					ch=fgetc(fp1);
					if(!(snin==11 || ch=='\n'))
					{
						seqname[snin]=ch;
						snin++;
					}

					if(ch=='\n')
					{
						ch=fgetc(fp1);
						break;
					}
				}//end ofwhile	

			}//end of if
		
			 if(ch!='\n')  
				fputc(ch,fp2);

			ch=fgetc(fp1);

			if(ch==EOF)
			{
				break;
			}
	 		
		}while(ch!='>');//end of while


		endpos=ftell(fp2);

		//int endpos=begpos;  
		
		
		seq=new char[endpos + 2];
		memset(seq,'\0',endpos + 2);

   
		begpos=0;
		fseek(fp2,begpos,SEEK_SET);

		ch1=fgetc(fp2);
		int seqin=0;

		while(ch1!=EOF)
		{
			seq[seqin]=ch1;
			seqin++;
			ch1=fgetc(fp2);
		}
		
		
		winobj->setWindowLen(winlen);
		winobj->setpcity(pcity);
		winobj->settherHold(thold*1.0);
		winobj->setIntronExon(exoncut,intcut);
	
	
		cout<<"Length:"<<strlen(seq)<<endl;
		
		winobj->doFt(seq,"log");

		 
       //put the coding regions in a separate file


	   for(it1=winobj->vecCRegion.begin();it1!=winobj->vecCRegion.end();it1++)
	   {

		  unsigned long actstartpos=it1->start-((winlen/2)-1);
		  unsigned long actendpos=it1->end +(winlen/2);
		  
		  if(actendpos>strlen(seq))
			actendpos=strlen(seq);

		  
		  seqs=new char[actendpos -actstartpos +2];
		  memset(seqs,'\0',actendpos -actstartpos +2);
		  
		  int j=0;
		  
		  for(int i=actstartpos -1;i<actendpos;i++,j++)
		     seqs[j]=seq[i];
		 
		  
		  
			int len=strlen(seqs);
		 		  
		  
			cout<<"Length:"<<len<<endl;
			winobj1->setWindowLen(strlen(seqs));
		  
			winobj1->setpcity(pcity);
			winobj1->settherHold(thold*1.0);
			winobj1->doFt(seqs,"log");
		  
			fprintf(fpc,">fourier cds%d %s [%d - %d] %d %0.4f\n",seqno,seqname,actstartpos,actendpos,(actendpos-actstartpos)+1,winobj1->x[0]);
			fprintf(fpc,"%s\n",seqs);

		
	   }//end of for
		
		fclose(fp2);
		remove("tmpseq.tmp");
		winobj->vecCRegion.clear();
		seqno++;
		
		delete []seq;
		seq = NULL;
		
   }//end of while
	  
	  
  fclose(fpc);

     
  }//end of main function

 
 int main(int argc, char* argv[])
 {
	
	ftcds  *Ch = new ftcds();

	int pcity;
	int excut,incut;
	
	cout<<argc<<endl;
	//cout<<"Usage:[filename] -w [length] -t [threshold]  -p [periodicity] -i [inputfile] -o [outputfile]"<<endl;
    		
	int winlen,thold;
	char *infile,*outfile;
	
	if(argc!=14)
	{
	   cout<<"Usage:[filename] -w [length] -t [threshold]  -p [periodicity] -ei [exoncut] [introncut] -i [inputfile] -o [outputfile]"<<endl;
	   exit(1);
	}

	for(int s=1;s<argc;s++)
    {
		
		if(strcmp(argv[s],"-w")==0)
		  winlen=atoi(argv[++s]);
   
		if(strcmp(argv[s],"-t")==0)
		 thold=atoi(argv[++s]);
   
		if(strcmp(argv[s],"-i")==0)
		{
			infile=new char[strlen(argv[s +1])+1];
			strcpy(infile,argv[++s]);
		}
   
		if(strcmp(argv[s],"-o")==0)
		{
			outfile=new char[strlen(argv[s +1])+1];
			strcpy(outfile,argv[++s]);
		}

		if(strcmp(argv[s],"-p")==0)
		   pcity=atoi(argv[++s]);

		if(strcmp(argv[s],"-ei")==0)
		{
		  excut=atoi(argv[s+1]);
		  incut=atoi(argv[s +2]);

		}

	}//end of for
    
	cout<<"lenwindow:"<<winlen<<endl;
	cout<<"Pcity:"<<pcity<<endl;
	cout<<"Exon intron"<<excut<<" "<<incut<<endl;
	cout<<"Thold:"<<thold<<endl;
	cout<<"Infile:"<<infile<<endl;
	cout<<"Outfile:"<<outfile<<endl;

	Ch->doGenscan(infile,outfile,winlen,thold,pcity,excut,incut);

	//Ch->doGenscan("c:\\anift\\amrna.txt","c:\\anift\\amrnaout.txt",300,4,3,50,50);
    
	
    delete Ch;

	return 0;
    	
	
 }//end main
